home *** CD-ROM | disk | FTP | other *** search
/ Mac Mania 5 / MacMania 5.toast / / Tools&Utilities / Plotfoil 3.2 / naca.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-18  |  6.3 KB  |  297 lines  |  [TEXT/MMCC]
  1. /*
  2.  * naca.c 
  3.  *
  4.  * Created by Alexandre Naaman (hoser@step.polymtl.ca)  30-08-1995
  5.  * Code clean-up, modifications to include more 5-digit sections: 
  6.  * Shamim Mohamed (shamim@synopsys.com) Sept. 8 95
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include <math.h>
  11. #include <stdlib.h>
  12. #include <string.h>
  13. #include <ctype.h>
  14.  
  15. #ifndef PI
  16. #define PI M_PI
  17. #endif
  18.  
  19. void
  20. usage(char *progname)
  21. {
  22.    fprintf(stderr, "Usage: %s [-n npoints] [-o output-file] <section-name>\n",
  23.        progname);
  24.  
  25. }
  26. enum series { four_digit, five_digit };
  27. enum surface { upper, lower };
  28.  
  29. struct airfoil {
  30.    char *name;
  31.    float maxor;
  32.    float posmax;
  33.    float thmax;
  34.    int k1;
  35.    enum series serie;
  36. };
  37.  
  38. float
  39. sqr(float x) {
  40.    return x*x; 
  41. }
  42.  
  43. float 
  44. cube(float x) {
  45.    return x*x*x; 
  46. }
  47.  
  48. float
  49. xcoor(float angl)
  50. {
  51.    return 0.5 + cos(angl)/2.0;
  52. }
  53.  
  54. float
  55. ytfunc(float x, float thmax)
  56. {
  57.    return thmax * (1.4845*sqrt(x) - 0.63*x - 1.758*sqr(x)+
  58.                1.4215*cube(x) - 0.5075*sqr(sqr(x)));
  59. }
  60.  
  61. void
  62. camber_four(float *yc, float *slope,
  63.       float x, float maxor, float posmax)
  64. {
  65.    if (x < posmax) {
  66.       *yc = maxor/sqr(posmax) * (2.0*posmax*x - sqr(x));
  67.       *slope = 2*(maxor/posmax) * (1.0-(x/posmax));
  68.    }
  69.    else {
  70.       *yc = maxor/sqr(1.0-posmax) * ( 1 - 2*posmax + 2*posmax*x - sqr(x));
  71.       *slope = 2*(maxor/sqr(1-posmax))*(posmax-x);
  72.    }
  73. }
  74.  
  75. void 
  76. camber_five(float *yc, float *slope,
  77.         float x, float maxor, float posmax, int k1)
  78. {
  79.    if (x < maxor) {
  80.       *yc = (1.0/6.0)*k1*(cube(x) - 3.0*maxor*sqr(x) + sqr(maxor)*(3.0-maxor)*x);
  81.       *slope = (1.0/6.0)*k1*( 3.0*sqr(x) - 6.0*maxor*x +
  82.                  sqr(maxor)*(3.0-maxor) );
  83.    }
  84.    else {
  85.       *yc = (1.0/6.0)*k1*cube(maxor)*(1.0-x);
  86.       *slope = (1.0/6.0)*k1*cube(maxor);
  87.    }
  88. }
  89.  
  90. void 
  91. out_point(float x, float yc, float yt, float slope, enum surface side,
  92.       FILE *fp)
  93. {
  94.    float xloc, yloc, h;
  95.  
  96. #ifdef Debug
  97.    fprintf(stderr, " >>> x: %.4g\tcamber: %.3g\tthickness: %.3g\n", x, yc, yt);
  98. #endif
  99.  
  100.    h = sqrt(1+sqr(slope));
  101.    if (side == upper) {
  102.       xloc = x - yt*slope/h;
  103.       yloc = yc + yt/h;
  104.    }
  105.    else {
  106.       xloc = x + yt*slope/h;
  107.       yloc = yc - yt/h;
  108.    }
  109.    /* Shove results into file */
  110.    fprintf(fp,"%5.5f %5.5f\n", xloc, yloc);
  111. }
  112.  
  113. #define dig(c) ((c)-'0')
  114.  
  115. /* ref. Abbott and Dozenthoff(?) */
  116. float cambers[4][5] = {
  117. { 1.1, 1.5, 1.8, 2.1, 2.3 },
  118. { 0.0, 2.3, 2.8, 3.1, 0.0 },
  119. { 0.0, 3.1, 3.7, 4.2, 0.0 },
  120. { 0.0, 4.6, 5.5, 6.2, 0.0 } };
  121.  
  122. void 
  123. get_params(struct airfoil *data, const char *name)
  124. {
  125.    data->name = name;
  126.    if(strlen(name) == 4) {
  127.       data->serie = four_digit;
  128.       data->maxor = dig(name[0])/100.0;
  129.       data->posmax = dig(name[1])/10.0;
  130.       data->thmax = atoi(name+2)/100.0;
  131.    }
  132.    else {
  133.       int d1, d2;
  134.       float max_camber;
  135.  
  136.       data->serie = five_digit;
  137.       data->thmax = atoi(name+3)/100.0;
  138.       d1 = dig(name[0]); d2 = dig(name[1]);
  139.       data->posmax = d2/20.0;
  140.       if(name[2] == '0') {
  141.      max_camber = cambers[d1-2][d2-1];
  142.      /* what the hell is this k1? And the the max camber points don't
  143.      // seem to correspond with the this code.... */
  144.      data->maxor = max_camber;
  145.      /*
  146.      switch (code) {
  147.      case 210: data->maxor=.058; data->posmax=.05; data->k1=361.4; break;
  148.      case 220: data->maxor=.126; data->posmax=.1;  data->k1=51.64; break;
  149.      case 230: data->maxor=.2025;data->posmax=.15; data->k1=15.957;break;
  150.      case 240: data->maxor=.290; data->posmax=.2;  data->k1=6.643; break;
  151.      case 250: data->maxor=.391; data->posmax=.25; data->k1=3.23;  break;
  152.      default: break;
  153.      }
  154.      */
  155.       }
  156.       /* else ... */
  157.    }
  158. }      
  159.  
  160. void 
  161. draw_surface(int nbrdiv, const struct airfoil *data, enum surface side,
  162.          FILE *fp)
  163. {
  164.    float angl;
  165.    float ainc = PI / nbrdiv;
  166.    int curnbr = nbrdiv; 
  167.    float x, yt, yc, slope;
  168.  
  169.    if (side == upper) {
  170.       angl = PI; ainc = - PI / nbrdiv;
  171.    }
  172.    else  {
  173.       angl = 0.0; ainc = PI / nbrdiv;
  174.    }
  175.  
  176.    for(; curnbr-- > 0; angl += ainc) {
  177.       x = xcoor(angl);
  178.       yt = ytfunc(x, data->thmax);
  179.       if (data->serie == four_digit)
  180.      camber_four(&yc, &slope, x, data->maxor, data->posmax);
  181.       else
  182.      camber_five(&yc, &slope, x, data->maxor, data->posmax, data->k1); 
  183.       out_point(x, yc, yt, slope, side, fp);
  184.    }
  185. }
  186.  
  187. /*
  188.  * Check to see if this is a valid name for an NACA section
  189.  */
  190. char *
  191. check_name(char *name)
  192. {
  193.    int len;
  194.  
  195.    /* Trim off a leading NACA or naca */
  196.    if(!strncmp(name, "NACA", 4) || !strncmp(name, "naca", 4))
  197.       name += 4;
  198.  
  199.    len = strlen(name);
  200.  
  201.    if(len < 4 || len > 5) {
  202.       return 0;
  203.    }
  204.    
  205.    if(len == 4) {
  206.       /* Need to all be numeric */
  207.       int i;
  208.       for(i = 0; i < 4; i++) 
  209.      if(!isdigit(name[i]))
  210.         return 0;
  211.    }
  212.    else {
  213.       /* Ok, first digit has to be 2, 3, 4, 6 */
  214.       if(!index("2346", name[0]))
  215.      return 0;
  216.  
  217.       /* Second digit has to be 12345, except if the first digit is [346]
  218.       // then the second can only be [234] */
  219.       if(name[1] < '1' || name[1] > '5')
  220.      return 0;
  221.       if(index("346", name[0]))
  222.      if(!index("234", name[1]))
  223.         return 0;
  224.  
  225.       /* Third digit has to be 0 or 1 */
  226.       if(name[2] < '0' || name[2] > '1')
  227.      return 0;
  228.  
  229.       /* We don't know how to do reflexed camber... */
  230.       if(name[2] == '1') {
  231.      fputs("Sorry, don't know how to do reflex camber (yet).\n", stderr);
  232.      return 0;
  233.       }
  234.  
  235.       /* last two need to be numeric */
  236.       if(!isdigit(name[3]) || !isdigit(name[4]))
  237.      return 0;
  238.    }
  239.    return name;
  240. }
  241.   
  242. int 
  243. generate_naca(char *name, int num, FILE *fp)
  244. {
  245.    struct airfoil data;
  246.  
  247.    get_params(&data, name);
  248.    
  249.    fprintf(fp, "NACA %s\n", name);
  250. /*   fprintf(fp,"0.0 0.0\n"); */
  251.  
  252. #ifdef Debug
  253.    fprintf(stderr, "%s:\n  camber: %g\n   thickness: %g\n   max_camber_pos: %g\n   series: %s\n",
  254.        data.name, data.maxor, data.thmax, data.posmax,
  255.        ((data.serie == four_digit)? "four digit" : "five digit"));
  256. #endif
  257.     
  258.    draw_surface(num, &data, upper, fp);
  259.    draw_surface(num, &data, lower, fp);
  260.    
  261.    fprintf(fp, "0.0 0.0\n");
  262.    fclose (fp);
  263.    
  264.    return 0;
  265. }
  266.  
  267. int 
  268. main(int argc, char *argv[])
  269. {
  270.    int argp = 1;
  271.    FILE *fp = stdout;
  272.    int num_points = 20;
  273.    char *name;
  274.    
  275.    if(argc <= 1) {
  276.       usage(argv[0]);
  277.       exit(1);
  278.    }
  279.  
  280.    /*
  281.     * process options
  282.     * right now, we don't have any but ew need to do -n and -o
  283.     */
  284.  
  285.    if(!(name = check_name(argv[argp]))) {
  286.       fprintf(stderr, "%s: \"%s\" is not a valid NACA section name.\n",
  287.           argv[0], argv[argp]);
  288.       exit(0);
  289.    }
  290.  
  291.    /* Open the output file for write */
  292.  
  293.    generate_naca(name, num_points, fp);
  294.    exit(0);
  295. }
  296.  
  297.